%计算在外磁场B(对应附加能量系数b)、最近邻作用耦合系数J=1条件下并考虑周期边界条件时的体系哈密顿量（以k为单位，相当于无量纲化）
%E=[-0.5*J*sum_{i,j=1}^n \sigma_{ij} * 最近邻 - b*sum_{i,j=1}^n \sigma_{ij}]/T
%其中，J前面的0.5是因为遍历时会多算一遍
function E = ising_energy(mat, n, b)
    J = 1;
    E = - b * sum(sum(mat));
    for i = 1 : n
        for j = 1 : n
            if(i == 1)
                if(j==1)
                    sumi = mat(1, 2) + mat(2, 1) + mat(1, n) + mat(n, 1);
                elseif(j==n)
                    sumi = mat(1, n-1) + mat(2, n) + mat(n, n) + mat(1, 1);
                else
                    sumi = mat(1, j-1) + mat(1, j+1) + mat(2, j) + mat(n, j);
                end
            elseif(i == n)
                if(j==1)
                    sumi = mat(n, 2) + mat(n-1, 1) + mat(1, 1) + mat(n, n);
                elseif(j==n)
                    sumi = mat(n, n-1) + mat(n-1, n) + mat(1, n) + mat(n, 1);
                else
                    sumi = mat(n, j-1) + mat(n, j+1) + mat(n-1, j) + mat(1, j);
                end
            else
                if(j==1)
                    sumi = mat(i-1, 1) + mat(i+1, 1) + mat(i, 2) + mat(i, n);
                elseif(j==n)
                    sumi = mat(i-1, n) + mat(i+1, n) + mat(i, n-1) + mat(i, 1);
                else
                    sumi = mat(i-1, j) + mat(i+1, j) + mat(i, j-1) + mat(i, j+1);
                end
            end
            E = E - 0.5 * J * mat(i,j) * sumi;
        end
    end
    return
end